home *** CD-ROM | disk | FTP | other *** search
/ Macwelt 1 / Macwelt DVD 1.toast / Web-Publishing / HTML-Editoren / Alpha ƒ / Mode Examples / S-Example.s < prev    next >
Encoding:
Text File  |  2000-10-30  |  4.4 KB  |  149 lines
  1. ## -*-S-*-
  2.  # ==========================================================================
  3.  # S-Example.s
  4.  # 
  5.  # Distributed as an example of Alpha's S mode.
  6.  # 
  7.  # S mode is available in the Statistical Modes package at
  8.  # 
  9.  # <ftp://ftp.ucsd.edu/pub/alpha/>
  10.  # ==========================================================================
  11.  ##
  12.  
  13. ## 
  14.  # ==========================================================================
  15.  #  FILE: "ps03.1.s"
  16.  #                                    created: 11/23/98 {09:30:54 pm} 
  17.  #                                last update: 05/29/00 {06:18:14 pm} 
  18.  #  Description: 
  19.  # 
  20.  #  S-Plus source code necessary to answer question 1 of problem set 3. 
  21.  #  Summary statistics can be obtained by running source file
  22.  #  ps03.1-summary.s
  23.  #
  24.  #  This file also uses the S-Plus functions "latab.reg1()" and
  25.  #  "print.reg.out()" that will print regression output in LaTeX format.
  26.  # 
  27.  #  Author: Craig Barton Upright
  28.  #  E-mail: <cupright@princeton.edu>
  29.  #    mail: Princeton University, Department of Sociology
  30.  #          Princeton, New Jersey 08544
  31.  #     www: <http://www.princeton.edu/~cupright>
  32.  #  
  33.  #  Copyright (c) 1998-2000 Craig Barton Upright
  34.  # 
  35.  #  All rights reserved.
  36.  # ==========================================================================
  37.  ##
  38.  
  39. ### Preliminaries
  40.  
  41. # reading in and transforming data, attaching dataframe data01
  42.  
  43. data03.1    <-  read.table("../data/ship.dat")
  44.  
  45. # correcting a misnamed factor value in data03.1$year
  46.  
  47. for (i in 1:length(data03.1$year)){
  48.     data03.1$year <- as.vector(data03.1$year)
  49.     if(data03.1$year[i] == "1965-70") (data03.1$year[i] <- c("1965-69"))
  50. }
  51.  
  52. data03.1$year    <-  factor(data03.1$year)
  53.  
  54. # creating a linear variable for "year," using midpoint of construction year
  55.  
  56. year.linear    <-  as.vector(data03.1$year)
  57.  
  58. for (i in 1:length(year.linear)){
  59.  
  60.     if(year.linear[i] == "1960-64") (year.linear[i] <- 1962)
  61.     if(year.linear[i] == "1965-69") (year.linear[i] <- 1967)    
  62.     if(year.linear[i] == "1970-74") (year.linear[i] <- 1972)
  63.     if(year.linear[i] == "1975-79") (year.linear[i] <- 1977)
  64. }
  65.  
  66. data03.1$year.linear   <-  as.numeric(year.linear); rm(year.linear)    
  67.  
  68. options(contrasts = c("contr.treatment", "contr.treatment") )
  69.  
  70. attach(data03.1)
  71.  
  72. out3.11 <-  glm(incidents ~ type1 + year + offset(log(months)),
  73.                 x = T, y = T, poisson)
  74. out3.12 <-  glm(incidents ~ type1 + year.linear + offset(log(months)), 
  75.                 x = T, y = T, poisson)
  76. out3.13    <-  glm(incidents ~ type1 + year.linear + year.linear^2 
  77.                 + offset(log(months)), poisson, x = T, y = T)
  78. out3.14 <-  glm(incidents ~ type1 * year.linear + offset(log(months)), 
  79.                 x = T, y = T, poisson)
  80.  
  81. ### Figures
  82.         
  83. # generating a postcript file with residual diagnostics of model out3.11
  84.  
  85. postscript(file = "../figures/ps03-QQPLOT.ps", horizontal = F, height = 4)
  86.  
  87.     par(mfrow = c(1,2))
  88.         plot(out3.11$y, resid(out3.11), 
  89.                 xlab = "observed values", ylab = "deviance residuals")
  90.         abline(h = 2, lty = 2)
  91.         qqnorm(resid(out3.11), ylab = "deviance residuals")
  92.         qqline(resid(out3.11), lty = 2)
  93.         dev.off()
  94.  
  95. detach()
  96.  
  97.  
  98. ### Summary Output
  99.  
  100. # creating a file with summary statistics, latex tables
  101.  
  102. source("/u/cupright/S-Plus_library-cbu/print.reg.out.s")
  103.  
  104. source("/u/cupright/S-Plus_library-cbu/latab.reg1.s")
  105.  
  106. d     <- date()
  107.  
  108. sink    ("ps03.1.txt")
  109.  
  110. cat("# ps03.1.txt \n")
  111. cat("# \n")
  112. cat("# Summary statistics of regression objects produced by ps03.1.s \n")
  113. cat("# source(\"ps03.1-summary.s\") \n")
  114. cat("# \n")
  115. cat("# Craig Barton Upright\n")
  116. cat("#", substring(d,9,10), substring(d,4,7), ",", substring(d,25,28), "\n\n")
  117.  
  118. print.reg.out(out3.11, name = "out3.11", 
  119.                 alt = "additive log-linear model")
  120. print.reg.out(out3.12, name = "out3.12", 
  121.                 alt = "treating year as linear effect")
  122. print.reg.out(out3.13, name = "out3.13", 
  123.                 alt = "treating year as linear with square effect")
  124. print.reg.out(out3.14, name = "out3.14", alt = "interaction model")
  125.  
  126. attach(data03.1)
  127.  
  128. print(anova(out3.14, test = "F"))
  129.  
  130. detach()
  131.  
  132. cat("\n\nDeviance residuals for model out3.11: \n\n")
  133. print(round(resid(out3.11),2)); cat("\n\n")
  134.  
  135. cat("sum of deviance residuals for model out3.11:", 
  136.         round(sum(resid(out3.11)^2),2), "\n")
  137. cat("residual deviance of model out3.11:         ", 
  138.         round(out3.11$dev, 2), "\n\n")
  139.  
  140. latab.reg1(out3.11, dev = T)
  141. latab.reg1(out3.13, dev = T)
  142.  
  143. sink()
  144. rm(out3.11, out3.12, out3.13, out3.14)
  145. rm(d) 
  146. cat("\nGenerated ascii file ps03.1.txt \n\n")
  147.  
  148.  
  149. # .